*reads cross-sectional rasters in ASCII: cropland, gridded countries, pfaf4 converted to grid, irrigation, urban extent 
* these are from (1) original sources (2) R processing of .tif files, or (3) shapefile conversion to raster then ASCII in ArcGIS
* also create grid from location of flare 


*Cropland acreage in half degrees in 2000 from UBC
*follows read in cropland in R and aggregate to half degree

log using read_other_rasters, replace
timer on 1

clear

** work in the data directory until the end of this file
cd "../data"
 
ras2dta, file(Cropland2000) genxcoord(x) genycoord(y)  replace
summarize Cropland2000, detail
rename Cropland2000 cropland
label var cropland "Cropland share in 2000"
keep x y cropland
sort y x 
save cropland2000, replace

/*
*Download files from
*https://sedac.ciesin.columbia.edu/downloads/data/gpw-v4/gpw-v4-national-identifier-grid-rev11/gpw-v4-national-identifier-grid-rev11_30_min_asc.zip 
*and unzip (can't do within Stata because login required)
* ASCII file renamed country_grid_30min 
* need package "kountry" to get country identifiers for later merges
*/

**# gridded country data from GPW
clear
ras2dta, file(country_grid_30min) genxcoord(x) genycoord(y)  replace 
drop if y<31 | y>300
replace y=y-30
replace country_grid_30min=. if country_grid_30min==32767
kountry country_grid_30min, from(iso3n) to(iso3c)
rename _ISO3C_ country_grid
label var country_grid "Country (ISO3c)"
keep y x country_grid
sort y x
save country_grid_30min, replace







**# Irrigation (from FAO Aquastat)
** extent is complete 90 N to 90 S 
** Data in 5 minutes, so each degree is 12 cells,
*** Dropping 15 degrees from top is 15*12=180 cells from top and obs below (90+60)*12=1800


copy "https://firebasestorage.googleapis.com/v0/b/fao-aquastat.appspot.com/o/GIS%2Fgmia_v5_aeigw_pct_aei_asc.zip?alt=media&token=ee5bd6d1-c8e2-44fd-a4cf-58d7563abbf6" irrig_gw.zip, replace
copy "https://firebasestorage.googleapis.com/v0/b/fao-aquastat.appspot.com/o/GIS%2Fgmia_v5_aeisw_pct_aei_asc.zip?alt=media&token=38b7da06-ddfa-4500-90df-2db65e27c9a7" irrig_sw.zip, replace

foreach type in gw sw {
unzipfile irrig_`type'.zip, replace

clear
ras2dta, file(gmia_v5_aei`type'_pct_aei) genxcoord(x) genycoord(y) dropmiss replace
sort y x

tempfile i_`type'
save `i_`type''
erase irrig_`type'.zip
erase gmia_v5_aei`type'_pct_aei.asc
erase gmia_v5_aei`type'_pct_aei.dta
}

copy "https://firebasestorage.googleapis.com/v0/b/fao-aquastat.appspot.com/o/GIS%2Fgmia_v5_aei_pct_asc.zip?alt=media&token=e448ce53-296f-4756-90c1-75c87f74e569" irrig_all.zip, replace
unzipfile irrig_all.zip, replace
clear
ras2dta, file(gmia_v5_aei_pct) genxcoord(x) genycoord(y) dropmiss replace
sort y x

merge 1:1 y x using `i_gw', keep(master match) nogen
merge 1:1 y x using `i_sw', keep(master match) nogen


drop if y<181 | y>1800
replace y=y-180



*GW and SW as share of irrigated area. Multiply to make it share of total land
gen igw=0 if gmia_v5_aei_pct==0
gen isw=0 if gmia_v5_aei_pct==0
replace igw=(gmia_v5_aeigw_pct_aei*gmia_v5_aei_pct)/10000 if gmia_v5_aei_pct>0 & !missing(gmia_v5_aei_pct)
replace isw=(gmia_v5_aeisw_pct_aei*gmia_v5_aei_pct)/10000 if gmia_v5_aei_pct>0 & !missing(gmia_v5_aei_pct)


* aggregate from 5 min (=.083 degree) to .5 degree
gen x_6=floor((x-1)/6) + 1
gen y_6=floor((y-1)/6) + 1 

collapse (mean) irrig=gmia_v5_aei_pct irrig_gw=igw irrig_sw=isw, by(x_6 y_6)

*share, not percent
replace irrig=irrig/100

label var irrig "Irrigation (share all land)"
label var irrig_gw "Groundwater irrigation (share all land)"
label var irrig_sw "Groundwater irrigation (share all land)"



rename x_6 x
rename y_6 y

sort x y 
compress
save irrigation, replace
erase irrig_all.zip
erase gmia_v5_aei_pct.asc
erase gmia_v5_aei_pct.dta



**# Pfaf4 grid

*Grid from Hydro1k shapefiles of level 4 subbasins 
clear
ras2dta, file(pfaf4_id) genxcoord(x) genycoord(y) dropmiss
rename pfaf4_id conpfaf4
sort x y
save pfaf4_grid, replace
erase pfaf4_id.dta


**# Admin 1 data in .05 deg grid
*Grid from GADM shapefiles of level 1 admin units 
clear
ras2dta, file(admin1) genxcoord(x) genycoord(y) replace 
sort x y
save admin1, replace


*Convert flares data into grid

import excel "https://eogdata.mines.edu/global_flare_data/VIIRS_Global_flaring_d.7_slope_0.0298_2012-2016_web.xlsx", sheet("flares_upstream") firstrow case(lower) clear

tempfile upstream
save `upstream'

import excel "https://eogdata.mines.edu/global_flare_data/VIIRS_Global_flaring_d.7_slope_0.0298_2012-2016_web.xlsx", sheet("flares_downstream_oil") firstrow case(lower) clear

tempfile do
save `do'

import excel "https://eogdata.mines.edu/global_flare_data/VIIRS_Global_flaring_d.7_slope_0.0298_2012-2016_web.xlsx", sheet("flares_downstream_gas") firstrow case(lower) clear

append using `upstream' `do'

gen int y=int((90-latitude)/.05)+1
gen int x=int((180+longitude)/.05)+1


sort y x

collapse (firstnm) type (mean) detection_frequency_2012 detection_frequency_2016, by(y x)

*drop first 15 degrees (=15*20 cells)
drop if y<301 | y>3100
replace y=y-300

label var type "Flare type"
label var detection_frequency_2012 "Flare detection freq (2012)"
label var detection_frequency_2016 "Flare detection freq (2016)"

label data "Flare frequency (cell only present if flare indicated, rest zeros)"

save flares, replace


**# Global Rural Urban Mapping
*Download files from "https://sedac.ciesin.columbia.edu/gpw/app/?code=ECU&file=grumpv1&data=urextent&type=ascii&resolut=30" and unzip
*takes a long time to read in
* file exceeds Dataverse size limits and is not provided
clear
cd gl_grumpv1_urextent_ascii_30
ras2dta, file(glurextents) genxcoord(x) genycoord(y) dropmiss replace
* extent is 85 degree N to -58 degrees S
* lights and droughts from 75 degrees N to -60 S so drop 10 degrees 
** each obs is 30 sec, so 60*2 obs per degree

use glurextents, clear
drop if y<1201 | y>18000
replace y=y-1200

* aggregate from 30 sec to .05 degree
gen x_6=floor((x-1)/6) + 1
gen y_6=floor((y-1)/6) + 1 

gen urban=(glurextents==2) if !missing(glurextents)

* characterize as urban if anywhere in .05 is urban 
collapse (max) urban (mean) urbanshare=urban, by(x_6 y_6)

label var urban "Urban (from GRUMP)"
label var urbanshare "Urban share"

rename x_6 x
rename y_6 y

sort x y 
compress

*back in data directory
cd ..
save urban, replace
erase "./gl_grumpv1_urextent_ascii_30/glurextents.dta"


cd "../code"


timer off 1
timer list 1

log close